/**
* @version:1.0
* @author kevinnan.org.cn  2018212114@mail.hfut.edu.cn
* @date 2020/8/1
* @Content:频率域滤波
*/
#include "frequency_filter.hpp"


frequencyfilter::frequencyFiltering::frequencyFiltering(cv::Mat image){
	this->inputImage = image.clone();
	fourierTransformImage = fourierTransform(inputImage);

	cv::Mat magImage = getMagnitudeImage(fourierTransformImage);

	magImage = changeCenter(magImage);

	//计算幅值,并转换到对数尺度
	//取对数的目的是使那些振幅较低的成分相对高振幅成分得以拉高，以便观察掩盖在低幅噪声中的周期信号。
	cv::log(magImage+1, magImage);

	//归一化处理，用0-1之间的浮点数将矩阵变换为可视的图像格式
    cv::normalize(magImage, magImage, 0, 1, CV_MINMAX);

    cv::imshow("originImage", inputImage);
    cv::imshow("fourierTransformResult", magImage);

    // cv::imwrite("../images/originImage.jpg", inputImage);
    // cv::imwrite("../images/fourierTransformResult.jpg", magImage);


    auto rImage = inverseFourierTransform(fourierTransformImage);

    //cv::imwrite("../images/inverseFourierTransformResult.jpg", rImage);
    cv::imshow("inverseFourierTransformResult", rImage);
}


cv::Mat frequencyfilter::frequencyFiltering::inverseFourierTransform(const cv::Mat& fourierImage){
	cv::Mat dealtImage = fourierImage.clone();
    cv::Mat iDFT[] = {cv::Mat::zeros(dealtImage.size(), CV_32FC1), 
    				  cv::Mat::zeros(dealtImage.size(), CV_32FC1)};

    cv::idft(dealtImage, dealtImage);
    cv::split(dealtImage, iDFT);
    //cv::log(iDFT[0]+1, iDFT[0]);
    //cout<<iDFT[0](cv::Rect(0,0,20,20))<<endl;
 //    normalize(iDFT[0], iDFT[0], 225, 0, CV_MINMAX);
 //    iDFT[0] = graytransform::gammaTransform(iDFT[0], 0.4);
	// iDFT[0].convertTo(iDFT[0], CV_32FC1);
    normalize(iDFT[0], iDFT[0], 1, 0, CV_MINMAX);
    return iDFT[0];
}


cv::Mat frequencyfilter::frequencyFiltering::fourierTransform(const cv::Mat& originImage){
	cv::Mat paddingImage;

	//得到适合傅里叶变换的最优尺寸
	int m = cv::getOptimalDFTSize(originImage.rows);
	int n = cv::getOptimalDFTSize(originImage.cols);

	//填充，上侧和左侧不填充
	cv::copyMakeBorder(originImage, paddingImage, 0, m-originImage.rows, 
					   0, n-originImage.cols, cv::BORDER_CONSTANT, cv::Scalar(0));

	//双通道Mat
	cv::Mat planes[] = {cv::Mat_<float>(paddingImage), cv::Mat::zeros(paddingImage.size(), CV_32FC1)};
	cv::Mat mergeImage;
	cv::merge(planes, 2, mergeImage);

	//进行傅里叶变换
	cv::dft(mergeImage, mergeImage, cv::DFT_COMPLEX_OUTPUT);
	
	return mergeImage;
}


cv::Mat frequencyfilter::frequencyFiltering::getMagnitudeImage(const cv::Mat& fourierImage){
	cv::Mat fourierImageCopy = fourierImage.clone();
	cv::Mat planes[] = {cv::Mat::zeros(fourierImageCopy.size(), CV_32FC1), 
						cv::Mat::zeros(fourierImageCopy.size(), CV_32FC1)};

	cv::Mat magImage = planes[0].clone();
	cv::split(fourierImageCopy, planes);
	cv::magnitude(planes[0], planes[1], magImage);

	//cv::log(magImage+1, magImage);

	//如果有奇数行或列，则转换为偶数
	magImage = magImage(cv::Rect(0, 0, magImage.cols-(magImage.cols%2), magImage.rows-(magImage.rows%2)));


	return magImage;
}


cv::Mat frequencyfilter::frequencyFiltering::getPhaseImage(const cv::Mat& fourierImage){
	cv::Mat fourierImageCopy = fourierImage.clone();
	cv::Mat planes[] = {cv::Mat::zeros(fourierImageCopy.size(), CV_32FC1), 
						cv::Mat::zeros(fourierImageCopy.size(), CV_32FC1)};

	cv::Mat phaseImage = planes[0].clone();
	cv::split(fourierImageCopy, planes);
	cv::phase(planes[0], planes[1], phaseImage);
	//imshow("phase", phaseImage);
	return phaseImage;
}


cv::Mat frequencyfilter::frequencyFiltering::changeCenter(const cv::Mat& magImage){
	//重新排布傅里叶变换后的图像，使得原点位于图像中心
	int centerX = magImage.cols / 2;
	int centerY = magImage.rows / 2;

	cv::Mat magImageCopy = magImage.clone();
	cv::Mat planes[] = {cv::Mat::zeros(magImageCopy.size(), CV_32FC1), 
						cv::Mat::zeros(magImageCopy.size(), CV_32FC1)};

	cv::Mat mat1(magImageCopy, cv::Rect(0, 0, centerX, centerY));				//左上
	cv::Mat mat2(magImageCopy, cv::Rect(0, centerY, centerX, centerY));			//右上
	cv::Mat mat3(magImageCopy, cv::Rect(centerX, 0, centerX, centerY));				//左下
	cv::Mat mat4(magImageCopy, cv::Rect(centerX, centerY, centerX, centerY));	//右下

	//互换左上和右下
	cv::Mat tempImage;
	mat1.copyTo(tempImage);
	mat4.copyTo(mat1);
	tempImage.copyTo(mat4);

	//互换左下和右上
	mat2.copyTo(tempImage);
	mat3.copyTo(mat2);
	tempImage.copyTo(mat3);

	return magImageCopy;
}


cv::Mat frequencyfilter::frequencyFiltering::idealLowPassFilter(float D0){
	//得到幅度谱
	cv::Mat magImage = getMagnitudeImage(this->fourierTransformImage);
	//将幅度谱中心化
	magImage = changeCenter(magImage);		
	//得到相位谱
	cv::Mat phaseImage = getPhaseImage(this->fourierTransformImage);

	//设置中心坐标
	int centerX = magImage.rows / 2;
	int centerY = magImage.cols / 2;

	for(int x = 0; x < magImage.rows; x++){
		for(int y = 0; y < magImage.cols; y++){
			double distance = sqrt(pow((x-centerX),2)+pow(y-centerY,2));
			if(distance <= D0){

			}else{
				//在D0范围之外，都为0，即过滤掉
				magImage.at<float>(x,y) = 0;
			}
		}
	}

	cv::Mat planes[] = {cv::Mat::zeros(magImage.size(), CV_32FC1), 
						cv::Mat::zeros(magImage.size(), CV_32FC1)};
	cv::Mat mergeImage = magImage.clone();
	cv::Mat changeImage = magImage.clone();

	//去中心化，显示变换后的图像
	magImage = changeCenter(magImage);

	//利用幅度谱和相位谱，变换得到傅里叶变换中的实数和复数分量，用于傅里叶逆变换
	cv::polarToCart(magImage, phaseImage, planes[0], planes[1]);

	cv::merge(planes, 2, mergeImage);

	//傅里叶逆变换，得到经理想低通滤波器变换后的图像
	changeImage = inverseFourierTransform(mergeImage);

    cv::imshow("idealLowPassFilter", changeImage);
    //cv::imwrite("../images/滤波结果.jpg", changeImage);

    return changeImage;
}


cv::Mat frequencyfilter::frequencyFiltering::butterworthLowPassFilter(float D0, int n){
	//得到幅度谱
	cv::Mat magImage = getMagnitudeImage(this->fourierTransformImage);
	//将幅度谱中心化
	magImage = changeCenter(magImage);		
	//得到相位谱
	cv::Mat phaseImage = getPhaseImage(this->fourierTransformImage);

	//设置中心坐标
	int centerX = magImage.rows / 2;
	int centerY = magImage.cols / 2;

	//BLPF滤波器
	cv::Mat BLPF = cv::Mat::zeros(magImage.size(), CV_32FC1);
	for(int x = 0; x < BLPF.rows; x++){
		for(int y = 0; y < BLPF.cols; y++){
			double distance = sqrt(pow((x-centerX),2)+pow(y-centerY,2));
			BLPF.at<float>(x,y) = 1 / (1+pow(distance/D0,2*n));
		}
	}

	cv::multiply(magImage, BLPF, magImage);

	cv::Mat planes[] = {cv::Mat::zeros(magImage.size(), CV_32FC1), 
						cv::Mat::zeros(magImage.size(), CV_32FC1)};
	cv::Mat mergeImage = magImage.clone();
	cv::Mat changeImage = magImage.clone();

	//去中心化，显示变换后的图像
	magImage = changeCenter(magImage);

	//利用幅度谱和相位谱，变换得到傅里叶变换中的实数和复数分量，用于傅里叶逆变换
	cv::polarToCart(magImage, phaseImage, planes[0], planes[1]);

	cv::merge(planes, 2, mergeImage);

	//傅里叶逆变换，得到经理想低通滤波器变换后的图像
	changeImage = inverseFourierTransform(mergeImage);

    cv::imshow("butterworthLowPassFilter", changeImage);

	return changeImage;
}


cv::Mat frequencyfilter::frequencyFiltering::guassLowPassFilter(float D0){
	//得到幅度谱
	cv::Mat magImage = getMagnitudeImage(this->fourierTransformImage);
	//将幅度谱中心化
	magImage = changeCenter(magImage);		
	//得到相位谱
	cv::Mat phaseImage = getPhaseImage(this->fourierTransformImage);

	//设置中心坐标
	int centerX = magImage.rows / 2;
	int centerY = magImage.cols / 2;

	//GLPF滤波器
	cv::Mat GLPF = cv::Mat::zeros(magImage.size(), CV_32FC1);
	for(int x = 0; x < GLPF.rows; x++){
		for(int y = 0; y < GLPF.cols; y++){
			double distance = sqrt(pow((x-centerX),2)+pow(y-centerY,2));
			GLPF.at<float>(x,y) = exp(-pow(distance, 2)/(2*pow(D0,2)));
		}
	}

	cv::multiply(magImage, GLPF, magImage);

	cv::Mat planes[] = {cv::Mat::zeros(magImage.size(), CV_32FC1), 
						cv::Mat::zeros(magImage.size(), CV_32FC1)};
	cv::Mat mergeImage = magImage.clone();
	cv::Mat changeImage = magImage.clone();

	//去中心化，显示变换后的图像
	magImage = changeCenter(magImage);

	//利用幅度谱和相位谱，变换得到傅里叶变换中的实数和复数分量，用于傅里叶逆变换
	cv::polarToCart(magImage, phaseImage, planes[0], planes[1]);

	cv::merge(planes, 2, mergeImage);

	//傅里叶逆变换，得到经理想低通滤波器变换后的图像
	changeImage = inverseFourierTransform(mergeImage);

    cv::imshow("guassLowPaddFilter", changeImage);

	return changeImage;
}


cv::Mat frequencyfilter::frequencyFiltering::idealHighPassFilter(float D0){
	//得到幅度谱
	cv::Mat magImage = getMagnitudeImage(this->fourierTransformImage);
	//将幅度谱中心化
	magImage = changeCenter(magImage);		
	//得到相位谱
	cv::Mat phaseImage = getPhaseImage(this->fourierTransformImage);

	//设置中心坐标
	int centerX = magImage.rows / 2;
	int centerY = magImage.cols / 2;

	for(int x = 0; x < magImage.rows; x++){
		for(int y = 0; y < magImage.cols; y++){
			double distance = sqrt(pow((x-centerX),2)+pow(y-centerY,2));
			if(distance <= D0){
				//在D0范围之内，都为0，即过滤掉
				magImage.at<float>(x,y) = 0;
			}else{

			}
		}
	}
	// cv::log(magImage+1, magImage);
	// cv::normalize(magImage, magImage, 0, 1, CV_MINMAX);
	// cv::imshow("fff", magImage);

	cv::Mat planes[] = {cv::Mat::zeros(magImage.size(), CV_32FC1), 
						cv::Mat::zeros(magImage.size(), CV_32FC1)};
	cv::Mat mergeImage = magImage.clone();
	cv::Mat changeImage = magImage.clone();

	//去中心化，显示变换后的图像
	magImage = changeCenter(magImage);

	//利用幅度谱和相位谱，变换得到傅里叶变换中的实数和复数分量，用于傅里叶逆变换
	cv::polarToCart(magImage, phaseImage, planes[0], planes[1]);

	cv::merge(planes, 2, mergeImage);

	//傅里叶逆变换，得到经理想低通滤波器变换后的图像
	changeImage = inverseFourierTransform(mergeImage);

	// changeImage.convertTo(changeImage, 0);
	// changeImage = changeImage(cv::Rect(0,0,inputImage.rows,inputImage.cols));
	// cout<<changeImage.type()<<endl;

    cv::imshow("idealHighPassFilter", changeImage);
    //cv::imwrite("../images/滤波结果.jpg", changeImage);

    return changeImage;
}


cv::Mat frequencyfilter::frequencyFiltering::butterworthHighPassFilter(float D0, int n){
	//得到幅度谱
	cv::Mat magImage = getMagnitudeImage(this->fourierTransformImage);
	//将幅度谱中心化
	magImage = changeCenter(magImage);		
	//得到相位谱
	cv::Mat phaseImage = getPhaseImage(this->fourierTransformImage);

	//设置中心坐标
	int centerX = magImage.rows / 2;
	int centerY = magImage.cols / 2;

	//BLPF滤波器
	cv::Mat BLPF = cv::Mat::zeros(magImage.size(), CV_32FC1);
	for(int x = 0; x < BLPF.rows; x++){
		for(int y = 0; y < BLPF.cols; y++){
			double distance = sqrt(pow((x-centerX),2)+pow(y-centerY,2));
			BLPF.at<float>(x,y) = 1 / (1+pow(D0/distance,2*n));
		}
	}

	cv::multiply(magImage, BLPF, magImage);

	cv::Mat planes[] = {cv::Mat::zeros(magImage.size(), CV_32FC1), 
						cv::Mat::zeros(magImage.size(), CV_32FC1)};
	cv::Mat mergeImage = magImage.clone();
	cv::Mat changeImage = magImage.clone();

	//去中心化，显示变换后的图像
	magImage = changeCenter(magImage);

	//利用幅度谱和相位谱，变换得到傅里叶变换中的实数和复数分量，用于傅里叶逆变换
	cv::polarToCart(magImage, phaseImage, planes[0], planes[1]);

	cv::merge(planes, 2, mergeImage);

	//傅里叶逆变换，得到经理想低通滤波器变换后的图像
	changeImage = inverseFourierTransform(mergeImage);

    cv::imshow("butterworthHighPassFilter", changeImage);

	return changeImage;
}


cv::Mat frequencyfilter::frequencyFiltering::guassHighPassFilter(float D0){
	//得到幅度谱
	cv::Mat magImage = getMagnitudeImage(this->fourierTransformImage);
	//将幅度谱中心化
	magImage = changeCenter(magImage);		
	//得到相位谱
	cv::Mat phaseImage = getPhaseImage(this->fourierTransformImage);

	//设置中心坐标
	int centerX = magImage.rows / 2;
	int centerY = magImage.cols / 2;

	//GLPF滤波器
	cv::Mat GLPF = cv::Mat::zeros(magImage.size(), CV_32FC1);
	for(int x = 0; x < GLPF.rows; x++){
		for(int y = 0; y < GLPF.cols; y++){
			double distance = sqrt(pow((x-centerX),2)+pow(y-centerY,2));
			GLPF.at<float>(x,y) = 1- exp(-pow(distance, 2)/(2*pow(D0,2)));
		}
	}

	cv::multiply(magImage, GLPF, magImage);

	cv::Mat planes[] = {cv::Mat::zeros(magImage.size(), CV_32FC1), 
						cv::Mat::zeros(magImage.size(), CV_32FC1)};
	cv::Mat mergeImage = magImage.clone();
	cv::Mat changeImage = magImage.clone();

	//去中心化，显示变换后的图像
	magImage = changeCenter(magImage);

	//利用幅度谱和相位谱，变换得到傅里叶变换中的实数和复数分量，用于傅里叶逆变换
	cv::polarToCart(magImage, phaseImage, planes[0], planes[1]);

	cv::merge(planes, 2, mergeImage);

	//傅里叶逆变换，得到经理想低通滤波器变换后的图像
	changeImage = inverseFourierTransform(mergeImage);

    cv::imshow("guassHighPaddFilter", changeImage/2);

	return changeImage;
}


cv::Mat frequencyfilter::frequencyFiltering::homomorphicFilter(float low, float high, float D0){
	cv::Mat originImage = inputImage.clone();
	originImage.convertTo(originImage, CV_32FC1);
	//对原图像做对数变换，将乘积分离为加法 
	cv::log(originImage+1, originImage);
	//做傅里叶变换
	cv::Mat fourierImage = fourierTransform(originImage);
	//得到幅度谱
	cv::Mat magImage = getMagnitudeImage(fourierImage);
	//将幅度谱中心化
	magImage = changeCenter(magImage);		
	//得到相位谱
	cv::Mat phaseImage = getPhaseImage(fourierImage);

	//设置中心坐标
	int centerX = magImage.rows / 2;
	int centerY = magImage.cols / 2;

	//HLPF滤波器
	cv::Mat HLPF = cv::Mat::zeros(magImage.size(), CV_32FC1);

	for(int x = 0; x < HLPF.rows; x++){
		for(int y = 0; y < HLPF.cols; y++){
			double distance = sqrt(pow((x-centerX),2)+pow(y-centerY,2));
			HLPF.at<float>(x,y) = (high-low)*(1- exp(-pow(distance, 2)/(2*pow(D0,2))))+low;
		}
	}

	cv::multiply(magImage, HLPF, magImage);

	cv::Mat planes[] = {cv::Mat::zeros(magImage.size(), CV_32FC1), 
						cv::Mat::zeros(magImage.size(), CV_32FC1)};
	cv::Mat mergeImage = magImage.clone();
	cv::Mat changeImage = magImage.clone();

	//去中心化，显示变换后的图像
	magImage = changeCenter(magImage);

	//利用幅度谱和相位谱，变换得到傅里叶变换中的实数和复数分量，用于傅里叶逆变换
	cv::polarToCart(magImage, phaseImage, planes[0], planes[1]);

	cv::merge(planes, 2, mergeImage);

	//傅里叶逆变换，得到经理想低通滤波器变换后的图像
	changeImage = inverseFourierTransform(mergeImage);

	cv::exp(changeImage, changeImage);
	//cv::normalize(changeImage, changeImage, 0, 1, CV_MINMAX);
	changeImage.convertTo(changeImage, 0);
	//inputImage.convertTo(inputImage, CV_32FC1);
	cout<<inputImage.size()<<"\t"<<changeImage.size()<<endl;
	changeImage = changeImage(cv::Rect(0, 0, inputImage.cols, inputImage.rows));

	cv::imshow("originImage_H", inputImage);
	cv::imshow("homomorphicFilter", (inputImage + changeImage)/2);
	
	return changeImage;
} 

